#ifndef PARTITIONER_DEF
#define PARTITIONER_DEF

#include <algorithm>
#include <queue>
#include <stack>
#include <map>
#include <random>
#include <unordered_set>

#include "api/graph.hpp"

class Basic{

public:
    void sequence(std::queue<vid_t>& stream, vid_t v_num)
    {
        for (vid_t v = 0; v < v_num; v++)
            stream.push(v);
        assert(v_num == stream.size());
    }

    void random_stream(std::queue<vid_t>& stream, vid_t v_num)
    {
        // 没有与参数列表匹配的shuffle实例
        std::vector<vid_t> arr(v_num);
        for (vid_t v = 0; v < v_num; v++)
            arr[v] = v;
        
        // 使用random_device生成随机数种子
        std::random_device rd;
        // 使用默认随机数引擎和种子
        std::default_random_engine rng(rd());
        
        // 使用shuffle函数将数组打乱
        std::shuffle(arr.begin(), arr.end(), rng); 
        for (vid_t v = 0; v < v_num; v++)
        {
            stream.push(arr[v]);
            // std::cout << arr[v] << " ";
        }
        // std::cout << std::endl;
        assert(v_num == stream.size());
    }

    void DFS_degree(EdgeContainer* graph_csr, std::queue<vid_t>& stream, vid_t v_num, std::vector<vid_t>& graph_degree)
    {
        std::vector<bool> is_visited(v_num, false);
        std::stack<vid_t> sta;
        vid_t max_degree = 0;
        vid_t max_d_i;

        for(vid_t v_i = 0; v_i < v_num; v_i++)
        {
            if(graph_degree[v_i] > max_degree && is_visited[v_i] == false){
                max_degree = graph_degree[v_i];
                max_d_i = v_i;
            }
        }
        sta.push(max_d_i);
        is_visited[max_d_i] = true;

        while(!sta.empty())
        {
            vid_t sta_top = sta.top();
            stream.push(sta_top);
            sta.pop();
            std::vector<vid_t> t;
            for(auto a = graph_csr->adj_lists[sta_top].begin; a < graph_csr->adj_lists[sta_top].end; a++) {
                if(is_visited[a->neighbor] == false){
                    t.push_back(a->neighbor);
                    is_visited[a->neighbor] = true;
                }
            }
            sort(t.begin(),t.end(),[&](vid_t a, vid_t b){
                return graph_degree[a] < graph_degree[b];
            });
            for(vid_t i =0; i< t.size();i++)
            {
                sta.push(t[i]);
            }
        }
        printf("stream size: %zu    %.2f%%\n",stream.size(),stream.size()/(float)v_num * 100);
        for(vid_t v = 0; v < v_num; v++){
            if(is_visited[v] == false){
                stream.push(v);
                is_visited[v] = true;
            }
        }
        assert(v_num == stream.size());
    }
};


class Naive : public Basic{

public:

    void naive(vid_t v_num, eid_t e_num, std::vector<std::vector<vid_t>>& partition_vertex_array, float beta,  std::vector<vid_t>& outPL)
    {
        logstream(LOG_DEBUG) << "Start to Naive!" << std::endl;
        int partition_num = partition_vertex_array.size();

        vid_t max_volumn = beta * e_num / partition_num; // e_num已经是无向图后的边数了，不需要*2
        logstream(LOG_INFO) << "max_volumn: " << max_volumn << std::endl;

        clock_t partition_start = clock();
        int p_i = 0;
        vid_t beg_degree = 0;
        for(vid_t v = 0; v < v_num; v++)
        {
            if(outPL[v + 1] - beg_degree <= max_volumn)
            {
                partition_vertex_array[p_i].push_back(v);
                continue;
            }
            p_i++;
            assert(p_i < partition_num);

            partition_vertex_array[p_i].push_back(v);
            beg_degree = outPL[v];  
        }
        clock_t partition_end = clock();
        logstream(LOG_DEBUG) << "Partition time is " << (double)(partition_end - partition_start) / CLOCKS_PER_SEC << "s" << std::endl;
    }
};

class MPGP : public Basic{

public:
    void mpgp(EdgeContainer* graph_csr, vid_t v_num, std::vector<std::vector<vid_t>>& partition_vertex_array, float beta, std::vector<vid_t>& graph_degree, std::string stream_mode)
    {
        logstream(LOG_DEBUG) << "Start to MPGP." << std::endl;
        int partition_num = partition_vertex_array.size();

        std::vector<float> partition_common_score(partition_num, 0);
        std::vector<vid_t> partition_volumn(partition_num, 0);
        std::queue<vid_t> q;
        vid_t current_v_num = 0;
        vid_t current_volumn = 0;
        
        if(stream_mode == "sequence"){
            sequence(q, v_num);
        }
        else if(stream_mode == "random"){
            random_stream(q, v_num);
        }
        else if(stream_mode == "DFS_degree"){
            DFS_degree(graph_csr, q, v_num, graph_degree);
            logstream(LOG_INFO) << "DEF_degree finished." << std::endl;
        }

        clock_t partition_start = clock();
        vid_t iter = 0;
        vid_t last_v = 0;
        while(!q.empty())
        {
            if(current_v_num - last_v > 1000){
                // printf("\rProgress: %.2f%%", current_v_num/(float)v_num * 100);
                logstream(LOG_DEBUG) << "Progress:" << current_v_num/(float)v_num * 100 << "%." << std::endl;
                last_v = current_v_num;
            }
            iter++;

            // float c = beta * current_v_num / partition_num;
            float c = beta * current_volumn / partition_num;

            vid_t v_i = q.front();
            q.pop();
            std::fill(partition_common_score.begin(), partition_common_score.end(), 0);

            // #pragma omp parallel for schedule(static)
            for(int p_i = 0; p_i < partition_num; p_i++)
            {
                std::vector<vid_t> v_intersection;
                int ps2_score;
                ps2_score = PS2(&graph_csr->adj_lists[v_i], partition_vertex_array[p_i], v_intersection, graph_degree, partition_volumn[p_i]);

                int ps1_score;
                ps1_score = PS1(v_intersection, v_i, graph_csr, partition_num);

                // float w = 1 - partition_vertex_array[p_i].size() / c;     
                float w = 1 - partition_volumn[p_i] / c;
                partition_common_score[p_i] = ((float)ps1_score + (float)ps2_score + 1) * w; // 为什么要+1?
            }

            int maxp = max_element(partition_common_score.begin(), partition_common_score.end()) - partition_common_score.begin();
            // logstream(LOG_DEBUG) << "v_i, maxp: " << v_i << ", " << maxp << std::endl;
            partition_vertex_array[maxp].insert(lower_bound(partition_vertex_array[maxp].begin(), partition_vertex_array[maxp].end(), v_i), v_i);
            partition_volumn[maxp] += graph_degree[v_i];
            current_v_num++;
            current_volumn += graph_degree[v_i];
        }
        clock_t partition_end = clock();
        logstream(LOG_DEBUG) << "Partition time is " << (double)(partition_end - partition_start) / CLOCKS_PER_SEC << "s" << std::endl;
        
        std::cout << "Vertices in each partition: ";
        for(int p_i = 0; p_i < partition_num; p_i++)
            std::cout << p_i << ": " << partition_vertex_array[p_i].size() << ", ";
        std::cout << std::endl;

        std::cout << "Edges in each partition: ";
        for(int p_i = 0; p_i < partition_num; p_i++)
            std::cout << p_i << ": " << partition_volumn[p_i] << ", ";
        std::cout << std::endl;

        assert(current_v_num == v_num);
    }

    int PS2(AdjList* v1, const std::vector<vid_t>& v2, std::vector<vid_t>& v_intersection, const std::vector<vid_t>& graph_degree, vid_t volumn)
    {
        int sz = v1->end - v1->begin;
        std::vector<vid_t> v1_v;
        
        for(int i = 0; i < sz; i++)
        {
            if(v1_v.empty() || (v1->begin + i)->neighbor != (v1->begin + i - 1)->neighbor)
                v1_v.push_back((v1->begin + i)->neighbor);
        }
        myIntersection(v1_v, v2, v_intersection);

        return v_intersection.size();
    }

    int PS1(const std::vector<vid_t>& vertex_array, vid_t v, EdgeContainer* cn_csr, int partition_num)
    {
        vid_t sum_cn = 0;
        // #pragma omp parallel for schedule(static)
        for(vid_t i = 0; i < vertex_array.size(); i++)
        {
            sum_cn += getCN_bs(v, vertex_array[i], cn_csr);
        }
        return sum_cn;
    }

    vid_t getCN_bs(vid_t src, vid_t dst, EdgeContainer* csr) 
    {
        AdjList* adj_list = &csr->adj_lists[src];
        AdjUnit target;
        target.neighbor = dst;
        auto edge = std::lower_bound(adj_list->begin, adj_list->end, target, [](const AdjUnit & a, const AdjUnit & b) { return a.neighbor < b.neighbor; });
        if (edge != adj_list->end && dst == edge->neighbor)
            return 1; // edge->data
        else
            return 0;
    }

    void myIntersection(std::vector<vid_t> v1, const std::vector<vid_t>& v2, std::vector<vid_t>& v_intersection) // v1 = v1_v, v2 = partition_vertex_array
    {
        int p1 = 0;
        int p2 = 0;
        int v1_sz = v1.size(); // |N(v)|
        int v2_sz = v2.size(); // |u∈P_i|
        const std::vector<vid_t>* long_v;
        const std::vector<vid_t>* short_v;
        if(v1_sz > v2_sz){
            long_v = &v1;
            short_v = &v2;
        }else{
            long_v = &v2;
            short_v = &v1;
        }
        int max_sz = std::max(v1_sz, v2_sz);
        int min_sz = std::min(v1_sz, v2_sz);
        
        if(v1.empty() || v2.empty()) return;
        if((*long_v)[max_sz - 1] < (*short_v)[0] || (*long_v)[0] > (*short_v)[min_sz - 1]) return;
        
        while(p1 < max_sz && p2 < min_sz){
            int offset = 1;
            int last_p = offset;
            if((*short_v)[p2] < (*long_v)[p1]){
                p2++;
                continue;
            }
            while((*long_v)[p1+offset-1] < (*short_v)[p2]){
                offset = offset * 2;
                last_p = p1 + offset < max_sz ? offset : max_sz - p1;
                if(p1 + offset >= max_sz) break;
            }
            if((*long_v)[max_sz-1] < (*short_v)[p2]){
                p2++;
                break;
            }
            auto iter = lower_bound(long_v->begin()+(p1+offset/2), long_v->begin()+p1+last_p, (*short_v)[p2]); // 如果在 区间找到了。用于在指定区域内查找不小于目标值的第一个元素
            int t = iter - long_v->begin();
            if(*iter==(*short_v)[p2]){
                v_intersection.push_back((*short_v)[p2]);
                p2++;
                p1=t++;
            }else{
                p2++;
                p1=(p1+offset/2);
            }
        }
    }

};

class AGPFlag : public Basic
{
public:
	double eps;
	uint seed;
	uint *H[2]; 
	double *finalReserve;
	uint finalreserve_count;
    uint *U[2]; 
	uint *candidate_set[2];
	uint candidate_count[2];
	double *residue[2];

    // std::vector<std::unordered_set<vid_t>> par_set;
    std::vector<std::vector<vid_t>> par_inflag;

    AGPFlag(vid_t v_num, int p_num)
    {
        logstream(LOG_DEBUG) << "Initial AGPFlag!" << std::endl;
        candidate_count[0] = 0;
        candidate_count[1] = 0;
        H[0] = new uint[v_num];
        H[1] = new uint[v_num];
        U[0] = new uint[v_num];
        U[1] = new uint[v_num];
        candidate_set[0] = new uint[v_num];
        candidate_set[1] = new uint[v_num];
        residue[0] = new double[v_num];
        residue[1] = new double[v_num];
        finalReserve = new double[v_num];
        finalreserve_count = 0;
        for (uint i = 0; i < v_num; i++)
        {
            residue[0][i] = 0;
            residue[1][i] = 0;
            finalReserve[i] = 0;
            H[0][i] = 0;
            H[1][i] = 0;
            U[0][i] = 0;
            U[1][i] = 0;
            candidate_set[0][i] = 0;
            candidate_set[1][i] = 0;
        }
        par_inflag.resize(p_num);
        for(vid_t i = 0; i < par_inflag.size(); i++){
            par_inflag[i].resize(v_num);
        }
        for(vid_t i = 0; i < par_inflag.size(); i++){
            for(vid_t j = 0; j < par_inflag[i].size(); j++){
                par_inflag[i][j] = false;
            }
            
        }
    }

    ~AGPFlag()
    {
        delete[] H[0];
        delete[] H[1];
        delete[] U[0];
        delete[] U[1];
        delete[] residue[0];
        delete[] residue[1];
        delete[] finalReserve;
        delete[] candidate_set[0];
        delete[] candidate_set[1];
    }

    void agp(EdgeContainer* graph_csr, vid_t v_num, std::vector<std::vector<vid_t>>& partition_vertex_array, float beta, std::vector<vid_t>& graph_degree, std::vector<vid_t>& outPL, std::vector<eid_t>& outEL, std::string stream_mode)
    {
        int partition_num = partition_vertex_array.size();

        std::vector<float> partition_common_score(partition_num, 0);
        std::vector<vid_t> partition_volumn(partition_num, 0);
        std::queue<vid_t> q;
        vid_t current_v_num = 0;
        vid_t current_volumn = 0;

        if(stream_mode == "sequence"){
            sequence(q, v_num);
        }
        else if(stream_mode == "random"){
            random_stream(q, v_num);
        }
        else if(stream_mode == "DFS_degree"){
            DFS_degree(graph_csr, q, v_num, graph_degree);
        }

        clock_t partition_start = clock();
        vid_t iter = 0;
        vid_t last_v = 0;
        float score;
        while(!q.empty())
        {
            if(current_v_num - last_v > 1000){
                logstream(LOG_DEBUG) << "Progress:" << current_v_num/(float)v_num * 100 << "%." << std::endl;
                last_v = current_v_num;
            }
            iter++;

            float c = beta * current_v_num / partition_num;
            // float c = beta * current_volumn / partition_num;
            logstream(LOG_INFO) << "current_volumn: " << current_volumn << ", c: " << c << std::endl;

            vid_t v_i = q.front();
            q.pop();
            std::fill(partition_common_score.begin(), partition_common_score.end(), 0);

            // #pragma omp parallel for schedule(static)
                for(int p_i = 0; p_i < partition_num; p_i++)
                {
                    if (partition_vertex_array[p_i].size() == 0) score = 1;
                    else score = compute_exp_length(graph_csr, v_i, p_i, graph_degree, outPL, outEL);
                    logstream(LOG_INFO) << "compute_exp_length in p" << p_i << ": " << score << std::endl;
                    
                    float w = 1 - partition_vertex_array[p_i].size() / c;
                    // float w = 1 - partition_volumn[p_i] / c;        
                    partition_common_score[p_i] = score * w;
                    logstream(LOG_INFO) << "partition_volumn[" << p_i << "]: " << partition_volumn[p_i] << ", partition_common_score[" << p_i << "]: " << partition_common_score[p_i] << std::endl;
                }
            int maxp = max_element(partition_common_score.begin(), partition_common_score.end()) - partition_common_score.begin();
            logstream(LOG_DEBUG) << "v_i, maxp: " << v_i << ", " << maxp << std::endl;
            partition_vertex_array[maxp].insert(lower_bound(partition_vertex_array[maxp].begin(), partition_vertex_array[maxp].end(), v_i), v_i); // 能否优化？
            par_inflag[maxp][v_i] = true;
            partition_volumn[maxp] += graph_degree[v_i];
            current_v_num++;
            current_volumn += graph_degree[v_i];
        }
        clock_t partition_end = clock();
        logstream(LOG_DEBUG) << "Partition time is " << (double)(partition_end - partition_start) / CLOCKS_PER_SEC << "s" << std::endl;

        std::cout << "Vertices in each partition: ";
        for(int p_i = 0; p_i < partition_num; p_i++)
            std::cout << p_i << ": " << partition_vertex_array[p_i].size() << ", ";
        std::cout << std::endl;

        std::cout << "Edges in each partition: ";
        for(int p_i = 0; p_i < partition_num; p_i++)
            std::cout << p_i << ": " << partition_volumn[p_i] << ", ";
        std::cout << std::endl;
        
        assert(current_v_num == v_num);
    }

    float compute_exp_length(EdgeContainer* graph_csr, vid_t u, int p, std::vector<vid_t>& degree, std::vector<vid_t>& outPL, std::vector<eid_t>& outEL)
    {
        std::vector<vid_t> inflag = par_inflag[p];

        // std::cout << "inflag: ";
        // for(vid_t i = 0; i < inflag.size(); i++)
        // {
        //     if(inflag[i] == true)
        //         std::cout << i << ", ";
        // }
        // std::cout << std::endl;

        for (uint j = 0; j < finalreserve_count; j++)
        {
            finalReserve[H[0][j]] = 0;
            H[1][H[0][j]] = 0;
        }
        
        finalreserve_count = 0;
        uint tempLevel = 0;

        float exp_length = 0;
        uint L = 4;
        std::vector<double> w_i(L+1, 0);
        w_i[L] = 1;
        double Y_i = 1;

        residue[0][u] = 1;
        candidate_set[0][0] = u;
        candidate_count[0] = 1;

        // bool flag = false;
        while (true)
        {
            uint tempLevelID = tempLevel % 2;
            uint newLevelID = (tempLevel + 1) % 2;
            uint candidateCnt = candidate_count[tempLevelID];
            if (candidateCnt == 0) break;
            candidate_count[tempLevelID] = 0;

            for (uint j = 0; j < candidateCnt; j++)
            {
                uint tempNode = candidate_set[tempLevelID][j];
                double tempR = residue[tempLevelID][tempNode];

                // logstream(LOG_INFO) << "From " << u << " to " << tempNode << "'s " << tempLevel << "-hop transition probability is：" << tempR << std::endl;
                exp_length += tempLevel * tempR;

                U[tempLevelID][tempNode] = 0;
                residue[tempLevelID][tempNode] = 0;
                if (H[1][tempNode] == 0)
                {
                    H[0][finalreserve_count++] = tempNode;
                    H[1][tempNode] = 1;
                }
                finalReserve[tempNode] += (w_i[tempLevel] / Y_i) * tempR; // 9: q^i(u) = q^i(u) + wi/Yi · r^i(u)
                if (Y_i < eps) continue;
                vid_t outSize = degree[tempNode];

                double incre = tempR * (1 - w_i[tempLevel] / Y_i) / outSize;
                if (tempR < eps * outSize) continue;
                AdjList* v1 = &graph_csr->adj_lists[tempNode];
                if (incre > eps)
                {
                    for (uint k = 0; k < outSize; k++)
                    {
                        vid_t newNode = (v1->begin + k)->neighbor; // newNode不一定在分区中

                        if(newNode == u || inflag[newNode] == true) // set.count(newNode) == 1
                        {
                            residue[newLevelID][newNode] += incre;  // 5: r^(i+1)(v) = r^(i+1)(v) + incre
                            if (U[newLevelID][newNode] == 0)
                            {                       
                                U[newLevelID][newNode] = 1;
                                candidate_set[newLevelID][candidate_count[newLevelID]++] = newNode;
                            }
                        }   
                    }
                }
                else
                {
                    double sampling_pr = incre / eps; // 6:
                    std::minstd_rand0 generator(seed);
                    std::binomial_distribution<uint> distribution(outSize, (double)sampling_pr);
                    uint expn = distribution(generator);
                    if (expn == outSize)
                    {
                        for (uint k = 0; k < outSize; k++)
                        {
                            vid_t newNode = (v1->begin + k)->neighbor;

                            if(newNode == u || inflag[newNode] == true) // set.find(newNode) != set.end()
                            {
                                residue[newLevelID][newNode] += incre;  // 5: r^(i+1)(v) = r^(i+1)(v) + incre
                                if (U[newLevelID][newNode] == 0)
                                {                       
                                    U[newLevelID][newNode] = 1;
                                    candidate_set[newLevelID][candidate_count[newLevelID]++] = newNode;
                                }
                            }
                        }
                    }
                    else
                    {
                        for (uint ri = 0; ri < expn; ri++)
                        {
                            uint tmpran = (uint)(rand() * (outSize - ri - 1)) + 1; // R.drand()的含义是什么？

                            uint idchangefar = outPL[tempNode] + ri + tmpran;
                            uint idchangenear = outPL[tempNode] + ri;
                            uint tmpchange = outEL[idchangefar];
                            outEL[idchangefar] = outEL[idchangenear];
                            outEL[idchangenear] = tmpchange;

                            vid_t newNode = (v1->begin + ri)->neighbor;

                            if(newNode == u || inflag[newNode] == true)
                            {
                                residue[newLevelID][newNode] += eps;  // 5: r^(i+1)(v) = r^(i+1)(v) + incre
                                if (U[newLevelID][newNode] == 0)
                                {                       
                                    U[newLevelID][newNode] = 1;
                                    candidate_set[newLevelID][candidate_count[newLevelID]++] = newNode;
                                }
                            }
                        }
                    }

                }
            }
            Y_i -= w_i[tempLevel];
            tempLevel++;
        }    
        return exp_length;
    }
};


















































// ------------------------------------------------------------------------------------------------------------------------
// ------------------------------------------------------------------------------------------------------------------------
class AGPMap : public Basic{

public:
	double eps;
	uint seed;
	uint *H[2]; 
	double *finalReserve;
	uint finalreserve_count;
    uint *U[2]; 
	uint *candidate_set[2];
	uint candidate_count[2];
	double *residue[2];

    std::vector<std::map<vid_t, vid_t>> par_v2id;
    std::vector<std::map<vid_t, vid_t>> par_id2v;

    AGPMap(vid_t v_num, int p_num)
    {
        logstream(LOG_DEBUG) << "Initial AGPMap!" << std::endl;
        candidate_count[0] = 0;
        candidate_count[1] = 0;
        H[0] = new uint[v_num];
        H[1] = new uint[v_num];
        U[0] = new uint[v_num];
        U[1] = new uint[v_num];
        candidate_set[0] = new uint[v_num];
        candidate_set[1] = new uint[v_num];
        residue[0] = new double[v_num];
        residue[1] = new double[v_num];
        finalReserve = new double[v_num];
        finalreserve_count = 0;
        for (uint i = 0; i < v_num; i++)
        {
            residue[0][i] = 0;
            residue[1][i] = 0;
            finalReserve[i] = 0;
            H[0][i] = 0;
            H[1][i] = 0;
            U[0][i] = 0;
            U[1][i] = 0;
            candidate_set[0][i] = 0;
            candidate_set[1][i] = 0;
        }
        par_v2id.resize(p_num);
        par_id2v.resize(p_num);
    }

    ~AGPMap()
    {
        delete[] H[0];
        delete[] H[1];
        delete[] U[0];
        delete[] U[1];
        delete[] residue[0];
        delete[] residue[1];
        delete[] finalReserve;
        delete[] candidate_set[0];
        delete[] candidate_set[1];
    }

    void agp(EdgeContainer* graph_csr, vid_t v_num, std::vector<std::vector<vid_t>>& partition_vertex_array, float beta, std::vector<vid_t>& graph_degree, std::vector<vid_t>& outPL, std::vector<eid_t>& outEL, std::string stream_mode)
    {
        int partition_num = partition_vertex_array.size();

        std::vector<float> partition_common_score(partition_num, 0);
        std::vector<vid_t> partition_volumn(partition_num, 0);
        std::queue<vid_t> q;
        vid_t current_v_num = 0;

        if(stream_mode == "sequence"){
            sequence(q, v_num);
        }
        else if(stream_mode == "random"){
            random_stream(q, v_num);
        }
        else if(stream_mode == "DFS_degree"){
            DFS_degree(graph_csr, q, v_num, graph_degree);
        }

        clock_t partition_start = clock();
        vid_t iter = 0;
        vid_t last_v = 0;
        float score;
        while(!q.empty())
        {
            if(current_v_num - last_v > 1000){
                logstream(LOG_DEBUG) << "Progress:" << current_v_num/(float)v_num * 100 << "%." << std::endl;
                last_v = current_v_num;
            }
            iter++;

            float c = beta * current_v_num / partition_num;
            vid_t v_i = q.front();
            q.pop();
            std::fill(partition_common_score.begin(), partition_common_score.end(), 0);

            // #pragma omp parallel for schedule(static)
                for(int p_i = 0; p_i < partition_num; p_i++)
                {
                    insertMap(p_i, v_i);
                    if (partition_vertex_array[p_i].size() == 0) score = 1;
                    else score = compute_exp_length(graph_csr, v_i, p_i, graph_degree, outPL, outEL);
                    // std::cout << "Vertex " << v_i << " in Partition " << p_i << "'s expected walk length is " << score << std::endl;
                    
                    float w = 1 - partition_vertex_array[p_i].size() / c;
                    partition_common_score[p_i] = score * w;
                    // logstream(LOG_INFO) << "partition_common_score[" << p_i << "]: " << partition_common_score[p_i] << std::endl;
                }

            int maxp = max_element(partition_common_score.begin(), partition_common_score.end()) - partition_common_score.begin();
            // logstream(LOG_DEBUG) << "v_i, maxp: " << v_i << ", " << maxp << std::endl;
            partition_vertex_array[maxp].insert(lower_bound(partition_vertex_array[maxp].begin(), partition_vertex_array[maxp].end(), v_i), v_i);
            for(int p_i = 0; p_i < partition_num; p_i++)
            {
                if(p_i != maxp)
                    deleteMap(p_i, v_i);
            }
            partition_volumn[maxp] += graph_degree[v_i];
            current_v_num++;
        }
        clock_t partition_end = clock();
        logstream(LOG_DEBUG) << "Partition time is " << (double)(partition_end - partition_start) / CLOCKS_PER_SEC << "s" << std::endl;
        assert(current_v_num == v_num);
    }

    float compute_exp_length(EdgeContainer* graph_csr, vid_t u, int p, std::vector<vid_t>& degree, std::vector<vid_t>& outPL, std::vector<eid_t>& outEL)
    {
        std::map<vid_t, vid_t> v2id = par_v2id[p];
        std::map<vid_t, vid_t> id2v = par_id2v[p];

        // std::cout << "v2id: ";
        // for(auto it:v2id)
        // {
        //     std::cout << it.first << ": " << it.second << ", ";
        // }
        // std::cout << std::endl;

        // std::cout << "id2v: ";
        // for(auto it:id2v)
        // {
        //     std::cout << it.first << ": " << it.second << ", ";
        // }
        // std::cout << std::endl;

        for (uint j = 0; j < finalreserve_count; j++)
        {
            finalReserve[H[0][j]] = 0;
            H[1][H[0][j]] = 0;
        }
        
        finalreserve_count = 0;
        uint tempLevel = 0;

        float exp_length = 0;
        uint L = 4;
        std::vector<double> w_i(L+1, 0);
        w_i[L] = 1;
        double Y_i = 1;
        residue[0][v2id[u]] = 1;
        candidate_set[0][0] = v2id[u];
        candidate_count[0] = 1;

        // bool flag = false;
        while (true)
        {
            uint tempLevelID = tempLevel % 2;
            uint newLevelID = (tempLevel + 1) % 2;
            uint candidateCnt = candidate_count[tempLevelID];
            if (candidateCnt == 0) break;
            candidate_count[tempLevelID] = 0;

            for (uint j = 0; j < candidateCnt; j++)
            {
                uint tempNodeID = candidate_set[tempLevelID][j];
                double tempR = residue[tempLevelID][tempNodeID];
                // logstream(LOG_INFO) << "From " << u << "(" << v2id[u] << ")" << " to " << id2v[tempNodeID] << "("  << tempNodeID << ")'s " << tempLevel << "-hop transition probability is：" << tempR << std::endl;
                exp_length += tempLevel * tempR;

                U[tempLevelID][tempNodeID] = 0;
                residue[tempLevelID][tempNodeID] = 0;
                if (H[1][tempNodeID] == 0)
                {
                    H[0][finalreserve_count++] = tempNodeID;
                    H[1][tempNodeID] = 1;
                }
                finalReserve[tempNodeID] += (w_i[tempLevel] / Y_i) * tempR; // 9: q^i(u) = q^i(u) + wi/Yi · r^i(u)
                if (Y_i < eps) continue;
                vid_t outSize = degree[id2v[tempNodeID]];

                double incre = tempR * (1 - w_i[tempLevel] / Y_i) / outSize;
                if (tempR < eps * outSize) continue;
                AdjList* v1 = &graph_csr->adj_lists[id2v[tempNodeID]];
                if (incre > eps)
                {
                    for (uint k = 0; k < outSize; k++)
                    {
                        vid_t newNode = (v1->begin + k)->neighbor; // newNode不一定在分区中

                        if(v2id.count(newNode) != 0)
                        {
                            vid_t newNodeId = v2id[newNode];
                            residue[newLevelID][newNodeId] += incre;  // 5: r^(i+1)(v) = r^(i+1)(v) + incre
                            if (U[newLevelID][newNodeId] == 0)
                            {                       
                                U[newLevelID][newNodeId] = 1;
                                candidate_set[newLevelID][candidate_count[newLevelID]++] = newNodeId;
                            }
                        }   
                    }
                }
                else
                {
                    double sampling_pr = incre / eps; // 6:
                    std::minstd_rand0 generator(seed);
                    std::binomial_distribution<uint> distribution(outSize, (double)sampling_pr);
                    uint expn = distribution(generator);
                    if (expn == outSize)
                    {
                        for (uint k = 0; k < outSize; k++)
                        {
                            vid_t newNode = (v1->begin + k)->neighbor;

                            if(v2id.count(newNode) != 0)
                            {
                                vid_t newNodeId = v2id[newNode];
                                residue[newLevelID][newNodeId] += incre;  // 5: r^(i+1)(v) = r^(i+1)(v) + incre
                                if (U[newLevelID][newNodeId] == 0)
                                {                       
                                    U[newLevelID][newNodeId] = 1;
                                    candidate_set[newLevelID][candidate_count[newLevelID]++] = newNodeId;
                                }
                            }
                        }
                    }
                    else
                    {
                        for (uint ri = 0; ri < expn; ri++)
                        {
                            uint tmpran = (uint)(rand() * (outSize - ri - 1)) + 1; // R.drand()的含义是什么？

                            uint idchangefar = outPL[id2v[tempNodeID]] + ri + tmpran;
                            uint idchangenear = outPL[id2v[tempNodeID]] + ri;
                            uint tmpchange = outEL[idchangefar];
                            outEL[idchangefar] = outEL[idchangenear];
                            outEL[idchangenear] = tmpchange;

                            vid_t newNode = (v1->begin + ri)->neighbor;

                            if(v2id.count(newNode) != 0)
                            {
                                vid_t newNodeId = v2id[newNode];
                                residue[newLevelID][newNodeId] += eps;  // 5: r^(i+1)(v) = r^(i+1)(v) + incre
                                if (U[newLevelID][newNodeId] == 0)
                                {                       
                                    U[newLevelID][newNodeId] = 1;
                                    candidate_set[newLevelID][candidate_count[newLevelID]++] = newNodeId;
                                }
                            }
                        }
                    }

                }
            }
            Y_i -= w_i[tempLevel];
            tempLevel++;
        }    
        return exp_length;
    }

    void testMap()
    {
        insertMap(0, 0);

        insertMap(0, 4); insertMap(1, 4); deleteMap(0, 4);

        insertMap(0, 3); insertMap(1, 3); deleteMap(1, 3);

        insertMap(0, 1); insertMap(1, 1); deleteMap(0, 1);

        insertMap(0, 2); insertMap(1, 2); deleteMap(1, 2);
    }

    void insertMap(int p, vid_t v)
    {
        // logstream(LOG_DEBUG) << "Insert vertex " << v << " to partition " << p << "." << std::endl; 
        std::map<vid_t, vid_t> v2id = par_v2id[p];
        std::map<vid_t, vid_t> id2v = par_id2v[p];

        assert(v2id.size() == id2v.size());

        if (v2id.size() == 0)
        {
            v2id[v] = 0;
            id2v[0] = v;
            par_v2id[p] = v2id;
            par_id2v[p] = id2v;
            return;
        }

        auto lastVert = id2v.rbegin();
        vid_t last_id = lastVert->first; // 最大的id
        last_id++;
        v2id[v] = last_id;
        id2v[last_id] = v;

        par_v2id[p] = v2id;
        par_id2v[p] = id2v;
    }

    void deleteMap(int p, vid_t v)
    {
        // logstream(LOG_DEBUG) << "Delete vertex " << v << " from partition " << p << "." << std::endl; 
        std::map<vid_t, vid_t> v2id = par_v2id[p];
        std::map<vid_t, vid_t> id2v = par_id2v[p];

        vid_t id;

        auto vit = v2id.find(v);
        assert(vit != v2id.end());
        // if(vit != v2id.end())
        {
            id = v2id[v];
            v2id.erase(vit);
        }
        auto idit = id2v.find(id);
        assert(idit != id2v.end());
        // if (idit != id2v.end())
        {
            id2v.erase(idit);
        }
        par_v2id[p] = v2id;
        par_id2v[p] = id2v;  
    }

};

// ------------------------------------------------------------------------------------------------------------------------
// ------------------------------------------------------------------------------------------------------------------------
class AGPVec : public Basic{

public:
	double eps;
	uint seed;
	uint *H[2]; 
	double *finalReserve;
	uint finalreserve_count;
    uint *U[2]; 
	uint *candidate_set[2];
	uint candidate_count[2];
	double *residue[2];

    // std::vector<std::vector<vid_t, vid_t>> par_v2id;
    std::vector<std::vector<vid_t>> par_id2v;

    AGPVec(vid_t v_num, int p_num)
    {
        logstream(LOG_DEBUG) << "Initial AGPVec!" << std::endl;
        candidate_count[0] = 0;
        candidate_count[1] = 0;
        H[0] = new uint[v_num];
        H[1] = new uint[v_num];
        U[0] = new uint[v_num];
        U[1] = new uint[v_num];
        candidate_set[0] = new uint[v_num];
        candidate_set[1] = new uint[v_num];
        residue[0] = new double[v_num];
        residue[1] = new double[v_num];
        finalReserve = new double[v_num];
        finalreserve_count = 0;
        for (uint i = 0; i < v_num; i++)
        {
            residue[0][i] = 0;
            residue[1][i] = 0;
            finalReserve[i] = 0;
            H[0][i] = 0;
            H[1][i] = 0;
            U[0][i] = 0;
            U[1][i] = 0;
            candidate_set[0][i] = 0;
            candidate_set[1][i] = 0;
        }
        // par_v2id.resize(p_num);
        par_id2v.resize(p_num);
    }

    ~AGPVec()
    {
        delete[] H[0];
        delete[] H[1];
        delete[] U[0];
        delete[] U[1];
        delete[] residue[0];
        delete[] residue[1];
        delete[] finalReserve;
        delete[] candidate_set[0];
        delete[] candidate_set[1];
    }

    void agp(EdgeContainer* graph_csr, vid_t v_num, std::vector<std::vector<vid_t>>& partition_vertex_array, float beta, std::vector<vid_t>& graph_degree, std::vector<vid_t>& outPL, std::vector<eid_t>& outEL, std::string stream_mode)
    {
        int partition_num = partition_vertex_array.size();

        std::vector<float> partition_common_score(partition_num, 0);
        std::vector<vid_t> partition_volumn(partition_num, 0);
        std::queue<vid_t> q;
        vid_t current_v_num = 0;

        if(stream_mode == "sequence"){
            sequence(q, v_num);
        }
        else if(stream_mode == "random"){
            random_stream(q, v_num);
        }
        else if(stream_mode == "DFS_degree"){
            DFS_degree(graph_csr, q, v_num, graph_degree);
        }

        clock_t partition_start = clock();
        vid_t iter = 0;
        vid_t last_v = 0;
        float score;
        while(!q.empty())
        {
            if(current_v_num - last_v > 1000){
                logstream(LOG_DEBUG) << "Progress:" << current_v_num/(float)v_num * 100 << "%." << std::endl;
                last_v = current_v_num;
            }
            iter++;

            float c = beta * current_v_num / partition_num;
            vid_t v_i = q.front();
            q.pop();
            std::fill(partition_common_score.begin(), partition_common_score.end(), 0);

            // #pragma omp parallel for schedule(static)
                for(int p_i = 0; p_i < partition_num; p_i++)
                {
                    insertVec(p_i, v_i);
                    if (partition_vertex_array[p_i].size() == 0) score = 1;
                    else score = compute_exp_length(graph_csr, v_i, p_i, graph_degree, outPL, outEL);
                    // logstream(LOG_INFO) << "Vertex " << v_i << " in Partition " << p_i << "'s expected walk length is " << score << std::endl;
                    
                    float w = 1 - partition_vertex_array[p_i].size() / c;
                    partition_common_score[p_i] = score * w;
                    // logstream(LOG_INFO) << "partition_common_score[" << p_i << "]: " << partition_common_score[p_i] << std::endl;
                }

            int maxp = max_element(partition_common_score.begin(), partition_common_score.end()) - partition_common_score.begin();
            // logstream(LOG_DEBUG) << "v_i, maxp: " << v_i << ", " << maxp << std::endl;
            partition_vertex_array[maxp].insert(lower_bound(partition_vertex_array[maxp].begin(), partition_vertex_array[maxp].end(), v_i), v_i);
            for(int p_i = 0; p_i < partition_num; p_i++)
            {
                if(p_i != maxp)
                    deleteVec(p_i, v_i);
            }
            partition_volumn[maxp] += graph_degree[v_i];
            current_v_num++;
        }
        clock_t partition_end = clock();
        logstream(LOG_DEBUG) << "Partition time is " << (double)(partition_end - partition_start) / CLOCKS_PER_SEC << "s" << std::endl;
        assert(current_v_num == v_num);
    }

    float compute_exp_length(EdgeContainer* graph_csr, vid_t u, int p, std::vector<vid_t>& degree, std::vector<vid_t>& outPL, std::vector<eid_t>& outEL)
    {
        std::vector<vid_t> id2v = par_id2v[p];

        // std::cout << "id2v: ";
        // for(vid_t i = 0; i < id2v.size(); i++)
        // {
        //     std::cout << i << ": " << id2v[i] << ", ";
        // }
        // std::cout << std::endl;

        for (uint j = 0; j < finalreserve_count; j++)
        {
            finalReserve[H[0][j]] = 0;
            H[1][H[0][j]] = 0;
        }
        
        finalreserve_count = 0;
        uint tempLevel = 0;

        float exp_length = 0;
        uint L = 4;
        std::vector<double> w_i(L+1, 0);
        w_i[L] = 1;
        double Y_i = 1;

        std::vector<vid_t>::iterator itr = std::find(id2v.begin(), id2v.end(), u);
        assert(itr != id2v.end());
        vid_t idx_u = std::distance(id2v.begin(), itr);

        residue[0][idx_u] = 1;
        candidate_set[0][0] = idx_u;
        candidate_count[0] = 1;

        // bool flag = false;
        while (true)
        {
            uint tempLevelID = tempLevel % 2;
            uint newLevelID = (tempLevel + 1) % 2;
            uint candidateCnt = candidate_count[tempLevelID];
            if (candidateCnt == 0) break;
            candidate_count[tempLevelID] = 0;

            for (uint j = 0; j < candidateCnt; j++)
            {
                uint tempNodeID = candidate_set[tempLevelID][j];
                double tempR = residue[tempLevelID][tempNodeID];

                // logstream(LOG_INFO) << "From " << u << "(" << idx_u << ")" << " to " << id2v[tempNodeID] << "("  << tempNodeID << ")'s " << tempLevel << "-hop transition probability is：" << tempR << std::endl;
                exp_length += tempLevel * tempR;

                U[tempLevelID][tempNodeID] = 0;
                residue[tempLevelID][tempNodeID] = 0;
                if (H[1][tempNodeID] == 0)
                {
                    H[0][finalreserve_count++] = tempNodeID;
                    H[1][tempNodeID] = 1;
                }
                finalReserve[tempNodeID] += (w_i[tempLevel] / Y_i) * tempR; // 9: q^i(u) = q^i(u) + wi/Yi · r^i(u)
                if (Y_i < eps) continue;
                vid_t outSize = degree[id2v[tempNodeID]];

                double incre = tempR * (1 - w_i[tempLevel] / Y_i) / outSize;
                if (tempR < eps * outSize) continue;
                AdjList* v1 = &graph_csr->adj_lists[id2v[tempNodeID]];
                if (incre > eps)
                {
                    for (uint k = 0; k < outSize; k++)
                    {
                        vid_t newNode = (v1->begin + k)->neighbor; // newNode不一定在分区中

                        itr = std::find(id2v.begin(), id2v.end(), newNode);

                        if(itr != id2v.end())
                        {
                            vid_t newNodeId = std::distance(id2v.begin(), itr);
                            residue[newLevelID][newNodeId] += incre;  // 5: r^(i+1)(v) = r^(i+1)(v) + incre
                            if (U[newLevelID][newNodeId] == 0)
                            {                       
                                U[newLevelID][newNodeId] = 1;
                                candidate_set[newLevelID][candidate_count[newLevelID]++] = newNodeId;
                            }
                        }   
                    }
                }
                else
                {
                    double sampling_pr = incre / eps; // 6:
                    std::minstd_rand0 generator(seed);
                    std::binomial_distribution<uint> distribution(outSize, (double)sampling_pr);
                    uint expn = distribution(generator);
                    if (expn == outSize)
                    {
                        for (uint k = 0; k < outSize; k++)
                        {
                            vid_t newNode = (v1->begin + k)->neighbor;
                            itr = std::find(id2v.begin(), id2v.end(), newNode);

                            if(itr != id2v.end())
                            {
                                vid_t newNodeId = std::distance(id2v.begin(), itr);
                                residue[newLevelID][newNodeId] += incre;  // 5: r^(i+1)(v) = r^(i+1)(v) + incre
                                if (U[newLevelID][newNodeId] == 0)
                                {                       
                                    U[newLevelID][newNodeId] = 1;
                                    candidate_set[newLevelID][candidate_count[newLevelID]++] = newNodeId;
                                }
                            }
                        }
                    }
                    else
                    {
                        for (uint ri = 0; ri < expn; ri++)
                        {
                            uint tmpran = (uint)(rand() * (outSize - ri - 1)) + 1; // R.drand()的含义是什么？

                            uint idchangefar = outPL[id2v[tempNodeID]] + ri + tmpran;
                            uint idchangenear = outPL[id2v[tempNodeID]] + ri;
                            uint tmpchange = outEL[idchangefar];
                            outEL[idchangefar] = outEL[idchangenear];
                            outEL[idchangenear] = tmpchange;

                            vid_t newNode = (v1->begin + ri)->neighbor;
                            itr = std::find(id2v.begin(), id2v.end(), newNode);

                            if(itr != id2v.end())
                            {
                                vid_t newNodeId = std::distance(id2v.begin(), itr);
                                residue[newLevelID][newNodeId] += eps;  // 5: r^(i+1)(v) = r^(i+1)(v) + incre
                                if (U[newLevelID][newNodeId] == 0)
                                {                       
                                    U[newLevelID][newNodeId] = 1;
                                    candidate_set[newLevelID][candidate_count[newLevelID]++] = newNodeId;
                                }
                            }
                        }
                    }

                }
            }
            Y_i -= w_i[tempLevel];
            tempLevel++;
        }    
        return exp_length;
    }

    void testVec()
    {
        insertVec(0, 0);

        insertVec(0, 4); insertVec(1, 4); deleteVec(0, 4);

        insertVec(0, 3); insertVec(1, 3); deleteVec(1, 3);

        insertVec(0, 1); insertVec(1, 1); deleteVec(0, 1);

        insertVec(0, 2); insertVec(1, 2); deleteVec(1, 2);
    }

    void insertVec(int p, vid_t v)
    {
        // logstream(LOG_DEBUG) << "Insert vertex " << v << " to partition " << p << "." << std::endl; 
        std::vector<vid_t> id2v = par_id2v[p];

        if (id2v.size() == 0)
        {
            id2v.push_back(v);
            par_id2v[p] = id2v;
            return;
        }

        id2v.push_back(v);
        par_id2v[p] = id2v;
    }

    void deleteVec(int p, vid_t v)
    {
        // logstream(LOG_DEBUG) << "Delete vertex " << v << " from partition " << p << "." << std::endl; 
        std::vector<vid_t> id2v = par_id2v[p];

        // 删除最后一个元素
        id2v.pop_back();
        par_id2v[p] = id2v;  
    }

};

// ------------------------------------------------------------------------------------------------------------------------
// ------------------------------------------------------------------------------------------------------------------------
class AGPSet : public Basic
{
public:
	double eps;
	uint seed;
	uint *H[2]; 
	double *finalReserve;
	uint finalreserve_count;
    uint *U[2]; 
	uint *candidate_set[2];
	uint candidate_count[2];
	double *residue[2];

    std::vector<std::unordered_set<vid_t>> par_set;
    // std::vector<std::vector<vid_t>> par_id2v;

    AGPSet(vid_t v_num, int p_num)
    {
        logstream(LOG_DEBUG) << "Initial AGPSet!" << std::endl;
        candidate_count[0] = 0;
        candidate_count[1] = 0;
        H[0] = new uint[v_num];
        H[1] = new uint[v_num];
        U[0] = new uint[v_num];
        U[1] = new uint[v_num];
        candidate_set[0] = new uint[v_num];
        candidate_set[1] = new uint[v_num];
        residue[0] = new double[v_num];
        residue[1] = new double[v_num];
        finalReserve = new double[v_num];
        finalreserve_count = 0;
        for (uint i = 0; i < v_num; i++)
        {
            residue[0][i] = 0;
            residue[1][i] = 0;
            finalReserve[i] = 0;
            H[0][i] = 0;
            H[1][i] = 0;
            U[0][i] = 0;
            U[1][i] = 0;
            candidate_set[0][i] = 0;
            candidate_set[1][i] = 0;
        }
        par_set.resize(p_num);
        // par_id2v.resize(p_num);
    }

    ~AGPSet()
    {
        delete[] H[0];
        delete[] H[1];
        delete[] U[0];
        delete[] U[1];
        delete[] residue[0];
        delete[] residue[1];
        delete[] finalReserve;
        delete[] candidate_set[0];
        delete[] candidate_set[1];
    }

    void agp(EdgeContainer* graph_csr, vid_t v_num, eid_t e_num, std::vector<std::vector<vid_t>>& partition_vertex_array, float beta, std::vector<vid_t>& graph_degree, std::vector<vid_t>& outPL, std::vector<eid_t>& outEL, std::string stream_mode)
    {
        int partition_num = partition_vertex_array.size();

        std::vector<float> partition_common_score(partition_num, 0);
        std::vector<vid_t> partition_volumn(partition_num, 0);
        std::queue<vid_t> q;
        vid_t current_v_num = 0;

        if(stream_mode == "sequence"){
            sequence(q, v_num);
        }
        else if(stream_mode == "random"){
            random_stream(q, v_num);
        }
        else if(stream_mode == "DFS_degree"){
            DFS_degree(graph_csr, q, v_num, graph_degree);
        }

        clock_t partition_start = clock();
        vid_t iter = 0;
        vid_t last_v = 0;
        float score;
        while(!q.empty())
        {
            if(current_v_num - last_v > 1000){
                logstream(LOG_DEBUG) << "Progress:" << current_v_num/(float)v_num * 100 << "%." << std::endl;
                last_v = current_v_num;
            }
            iter++;

            // float c = beta * current_v_num / partition_num;
            vid_t v_i = q.front();
            q.pop();
            std::fill(partition_common_score.begin(), partition_common_score.end(), 0);

            // #pragma omp parallel for schedule(static)
                for(int p_i = 0; p_i < partition_num; p_i++)
                {
                    if (partition_vertex_array[p_i].size() == 0) score = 1;
                    else score = compute_exp_length(graph_csr, v_i, p_i, graph_degree, outPL, outEL);
                    // logstream(LOG_INFO) << "Vertex " << v_i << " in Partition " << p_i << "'s expected walk length is " << score << std::endl;
                    
                    float c = beta * partition_volumn[p_i] / e_num;
                    float w = 1 - partition_vertex_array[p_i].size() / c;
                    partition_common_score[p_i] = score * w;
                }
            int maxp = max_element(partition_common_score.begin(), partition_common_score.end()) - partition_common_score.begin();
            // logstream(LOG_DEBUG) << "v_i, maxp: " << v_i << ", " << maxp << std::endl;
            partition_vertex_array[maxp].insert(lower_bound(partition_vertex_array[maxp].begin(), partition_vertex_array[maxp].end(), v_i), v_i); // 能否优化？
            par_set[maxp].insert(v_i);
            partition_volumn[maxp] += graph_degree[v_i];
            current_v_num++;
        }
        clock_t partition_end = clock();
        logstream(LOG_DEBUG) << "Partition time is " << (double)(partition_end - partition_start) / CLOCKS_PER_SEC << "s" << std::endl;
        assert(current_v_num == v_num);
    }

    float compute_exp_length(EdgeContainer* graph_csr, vid_t u, int p, std::vector<vid_t>& degree, std::vector<vid_t>& outPL, std::vector<eid_t>& outEL)
    {
        std::unordered_set<vid_t> set = par_set[p];

        // std::cout << "set: ";
        // for(auto it:set)
        // {
        //     std::cout << it << ", ";
        // }
        // std::cout << std::endl;

        for (uint j = 0; j < finalreserve_count; j++)
        {
            finalReserve[H[0][j]] = 0;
            H[1][H[0][j]] = 0;
        }
        
        finalreserve_count = 0;
        uint tempLevel = 0;

        float exp_length = 0;
        uint L = 4;
        std::vector<double> w_i(L+1, 0);
        w_i[L] = 1;
        double Y_i = 1;

        residue[0][u] = 1;
        candidate_set[0][0] = u;
        candidate_count[0] = 1;

        // bool flag = false;
        while (true)
        {
            uint tempLevelID = tempLevel % 2;
            uint newLevelID = (tempLevel + 1) % 2;
            uint candidateCnt = candidate_count[tempLevelID];
            if (candidateCnt == 0) break;
            candidate_count[tempLevelID] = 0;

            for (uint j = 0; j < candidateCnt; j++)
            {
                uint tempNode = candidate_set[tempLevelID][j];
                double tempR = residue[tempLevelID][tempNode];

                // logstream(LOG_INFO) << "From " << u << " to " << tempNode << "'s " << tempLevel << "-hop transition probability is：" << tempR << std::endl;
                exp_length += tempLevel * tempR;

                U[tempLevelID][tempNode] = 0;
                residue[tempLevelID][tempNode] = 0;
                if (H[1][tempNode] == 0)
                {
                    H[0][finalreserve_count++] = tempNode;
                    H[1][tempNode] = 1;
                }
                finalReserve[tempNode] += (w_i[tempLevel] / Y_i) * tempR; // 9: q^i(u) = q^i(u) + wi/Yi · r^i(u)
                if (Y_i < eps) continue;
                vid_t outSize = degree[tempNode];

                double incre = tempR * (1 - w_i[tempLevel] / Y_i) / outSize;
                if (tempR < eps * outSize) continue;
                AdjList* v1 = &graph_csr->adj_lists[tempNode];
                if (incre > eps)
                {
                    for (uint k = 0; k < outSize; k++)
                    {
                        vid_t newNode = (v1->begin + k)->neighbor; // newNode不一定在分区中

                        if(newNode == u || set.count(newNode) == 1) // set.count(newNode) == 1
                        {
                            residue[newLevelID][newNode] += incre;  // 5: r^(i+1)(v) = r^(i+1)(v) + incre
                            if (U[newLevelID][newNode] == 0)
                            {                       
                                U[newLevelID][newNode] = 1;
                                candidate_set[newLevelID][candidate_count[newLevelID]++] = newNode;
                            }
                        }   
                    }
                }
                else
                {
                    double sampling_pr = incre / eps; // 6:
                    std::minstd_rand0 generator(seed);
                    std::binomial_distribution<uint> distribution(outSize, (double)sampling_pr);
                    uint expn = distribution(generator);
                    if (expn == outSize)
                    {
                        for (uint k = 0; k < outSize; k++)
                        {
                            vid_t newNode = (v1->begin + k)->neighbor;

                            if(newNode == u || set.count(newNode) == 1) // set.find(newNode) != set.end()
                            {
                                residue[newLevelID][newNode] += incre;  // 5: r^(i+1)(v) = r^(i+1)(v) + incre
                                if (U[newLevelID][newNode] == 0)
                                {                       
                                    U[newLevelID][newNode] = 1;
                                    candidate_set[newLevelID][candidate_count[newLevelID]++] = newNode;
                                }
                            }
                        }
                    }
                    else
                    {
                        for (uint ri = 0; ri < expn; ri++)
                        {
                            uint tmpran = (uint)(rand() * (outSize - ri - 1)) + 1; // R.drand()的含义是什么？

                            uint idchangefar = outPL[tempNode] + ri + tmpran;
                            uint idchangenear = outPL[tempNode] + ri;
                            uint tmpchange = outEL[idchangefar];
                            outEL[idchangefar] = outEL[idchangenear];
                            outEL[idchangenear] = tmpchange;

                            vid_t newNode = (v1->begin + ri)->neighbor;

                            if(newNode == u || set.count(newNode) == 1)
                            {
                                residue[newLevelID][newNode] += eps;  // 5: r^(i+1)(v) = r^(i+1)(v) + incre
                                if (U[newLevelID][newNode] == 0)
                                {                       
                                    U[newLevelID][newNode] = 1;
                                    candidate_set[newLevelID][candidate_count[newLevelID]++] = newNode;
                                }
                            }
                        }
                    }

                }
            }
            Y_i -= w_i[tempLevel];
            tempLevel++;
        }    
        return exp_length;
    }
};

#endif
